clc;
clear 
close all
addpath('..\')

pwd
X = xlsread('风电场预测.xlsx');

n_in = 8;  % 输入
n_out = 1 ; 
or_dim = size(X,2) ;       % 记录特征数据维度
num_samples = 3000;  
scroll_window = 1;  
[res] = data_collation(X, n_in, n_out, or_dim, scroll_window, num_samples);


% 训练集和测试集划分

num_size = 0.8;                              
num_train_s = round(num_size * num_samples); 

%% 归一化
P_train = res(1: num_train_s,1);
P_train = reshape(cell2mat(P_train)',n_in*or_dim,num_train_s);
T_train = res(1: num_train_s,2);
T_train = cell2mat(T_train)';

P_test = res(num_train_s+1: end,1);
P_test = reshape(cell2mat(P_test)',n_in*or_dim,num_samples-num_train_s);
T_test = res(num_train_s+1: end,2);
T_test = cell2mat(T_test)';


%  数据归一化
[p_train, ps_input] = mapminmax(P_train, 0, 1);
p_test = mapminmax('apply', P_test, ps_input);

[t_train, ps_output] = mapminmax(T_train, 0, 1);
t_test = mapminmax('apply', T_test, ps_output);

vp_train = reshape(p_train,n_in,or_dim,num_train_s);
vp_test = reshape(p_test,n_in,or_dim,num_samples-num_train_s);

vt_train = t_train;
vt_test = t_test;

f_ = [size(vp_train,1) size(vp_train,2)];
outdim = n_out;

%%  创建网络
% 创建输入层
%  网络参数设置
filterSize = 5;  
dropoutFactor = 0.006;
numBlocks = 2; 
numFilters = 30;    
layer = sequenceInputLayer(f_, Normalization = "rescale-symmetric", Name = "input");
 
% 创建网络图
lgraph = layerGraph(layer);
outputName = layer.Name;
 
% 建立网络结构 -- 残差块
for i = 1 : numBlocks
    % 膨胀因子
    dilationFactor = 2^(i-1);
 
    % 创建TCN正向支路
    layers = [
        convolution1dLayer(filterSize, numFilters, DilationFactor = dilationFactor, Padding = "causal", Name="conv1_" + i)  % 一维卷积层 
        layerNormalizationLayer                                                                                             % 层归一化
        spatialDropoutLayer(dropoutFactor)                                                                                  % 空间丢弃层
        convolution1dLayer(filterSize, numFilters, DilationFactor = dilationFactor, Padding = "causal")                     % 一维卷积层  
        layerNormalizationLayer                                                                                             % 层归一化
        reluLayer                                                                                                           % 激活层
        spatialDropoutLayer(dropoutFactor)                                                                                  % 空间丢弃层
        additionLayer(4, Name = "add_" + i)
    ];
 
    % 添加残差块到网络
    lgraph = addLayers(lgraph, layers);
 
    % 连接卷积层到残差块
    lgraph = connectLayers(lgraph, outputName, "conv1_" + i);
 
    % 创建 TCN反向支路flip网络结构
    Fliplayers = [
        FlipLayer("flip_" + i)                                                                                               % 反向翻转
        convolution1dLayer(1, numFilters, Name = "convSkip_"+i);                                                             % 反向残差连接
        convolution1dLayer(filterSize, numFilters, DilationFactor = dilationFactor, Padding = "causal", Name="conv2_" + i)   % 一维卷积层
        layerNormalizationLayer                                                                                              % 层归一化
        spatialDropoutLayer(dropoutFactor)                                                                                   % 空间丢弃层
        convolution1dLayer(filterSize, numFilters, DilationFactor = dilationFactor, Padding = "causal")                      % 一维卷积层
        layerNormalizationLayer                                                                                              % 层归一化
        reluLayer                                                                                                            % 激活层
        spatialDropoutLayer(dropoutFactor, Name="drop" + i)                                                                  % 空间丢弃层
    ];
 
    % 添加 flip 网络结构到网络
    lgraph = addLayers(lgraph, Fliplayers);
 
    % 连接 flip 卷积层到残差块
    lgraph = connectLayers(lgraph, outputName, "flip_" + i);
    lgraph = connectLayers(lgraph, "drop" + i, "add_" + i + "/in3");
    lgraph = connectLayers(lgraph, "convSkip_"+i, "add_" + i + "/in4");
    % 残差连接 -- 首层
    if i == 1
        % 建立残差卷积层
        % Include convolution in first skip connection.
        layer = convolution1dLayer(1,numFilters,Name="convSkip");

        lgraph = addLayers(lgraph,layer);
        lgraph = connectLayers(lgraph,outputName,"convSkip");
        lgraph = connectLayers(lgraph,"convSkip","add_" + i + "/in2");
    else
        lgraph = connectLayers(lgraph,outputName,"add_" + i + "/in2");
    end
    
    % Update layer output name.
    outputName = "add_" + i;
end

tempLayers = flattenLayer("Name","flatten");
lgraph = addLayers(lgraph,tempLayers);


tempLayers = [
    bilstmLayer(50,'Name','lstmlayer')
    fullyConnectedLayer(outdim,"Name","fc")
    regressionLayer("Name","regressionoutput")];
lgraph = addLayers(lgraph,tempLayers);

lgraph = connectLayers(lgraph,outputName,"flatten");
lgraph = connectLayers(lgraph,"flatten","lstmlayer");


%%  训练参数设置

options = trainingOptions('adam', ...                 % 优化算法Adam  
    'MaxEpochs', 300, ...                            % 最大训练次数
    'GradientThreshold', 1, ...                       % 梯度阈值
    'InitialLearnRate', 0.01, ...         % 初始学习率
    'LearnRateSchedule', 'piecewise', ...             % 学习率调整
    'LearnRateDropPeriod', 120, ...                   % 训练120次后开始调整学习率
    'LearnRateDropFactor',0.1, ...                    % 学习率调整因子
    'L2Regularization', 0.001, ...         % 正则化参数
    'ExecutionEnvironment', 'gpu',...                 % 训练环境
    'Verbose', 1, ...                                 % 关闭优化过程
    'Plots', 'none');                    % 画出曲线

% %%  训练模型
% net = trainNetwork(vp_train,vt_train,lgraph,options);
% %%  查看网络结构
% figure
% plot(lgraph)
% title("BiTCN-biLSTM网络模型")
% set(gcf,'color','w')
% %analyzeNetwork(net);% 查看网络结构
% %  预测
% t_sim1 = predict(net, vp_train); 
% t_sim2 = predict(net, vp_test); 
% 
% %  数据反归一化
% T_sim1 = mapminmax('reverse', t_sim1, ps_output);
% T_sim2 = mapminmax('reverse', t_sim2, ps_output);


%% Adaboost增强学习部分
%  权重初始化
tic

D = ones(1, num_train_s) / num_train_s;
M = num_train_s;
N = num_samples-num_train_s;
%  参数设置
K = 5;                       % 弱回归器个数


%  弱回归器回归
for i = 1 : K
            i
    %  创建模型
    clear net

    net = trainNetwork(vp_train, vt_train, lgraph,options);

    result1 = predict(net, vp_train); 
    result2 =  predict(net, vp_test); 
     % 数据格式转换
    E_sim1 = result1;% cell2mat将cell元胞数组转换为普通数组
    E_sim2 = result2;

    %  仿真预测
    t_sim1(i, :) = E_sim1';
    t_sim2(i, :) = E_sim2';

    %  预测误差
    Error(i, :) = t_sim1(i, :) - t_train;

    %  调整D值
    weight(i) = 0;
    for j = 1 : M
        if abs(Error(i, j)) > 0.02
            weight(i) = weight(i) + D(i, j);
            D(i + 1, j) = D(i, j) * 1.1;
        else
            D(i + 1, j) = D(i, j);
        end
    end

    %  弱回归器i权重
    weight(i) = 0.5 / exp(abs(weight(i)));

    %  D值归一化
    D(i + 1, :) = D(i + 1, :) / sum(D(i + 1, :));

end

%  强预测器预测
weight = weight / sum(weight);

%  强回归器输出结果
T_sim1 = zeros(1, M);
T_sim2 = zeros(1, N);

for i = 1 : K
    output1 = (weight(i) * t_sim1(i, :));
    output2 = (weight(i) * t_sim2(i, :));

    T_sim1 = output1 + T_sim1;
    T_sim2 = output2 + T_sim2;
end
toc
%  数据反归一化
T_sim1 = mapminmax('reverse', T_sim1, ps_output);
T_sim2 = mapminmax('reverse', T_sim2, ps_output);
T_sim1 = double(T_sim1);
T_sim2 = double(T_sim2);


% 指标计算
disp('…………测试集误差指标…………')
[mae2,rmse2,mape2,error2]=calc_error(T_test,T_sim2);
fprintf('\n')


figure('Position',[200,300,600,200])
plot(T_test);
hold on
plot(T_sim2)
legend('真实值','预测值')
title('预测集预测效果对比')
xlabel('样本点')
ylabel('发电功率')


